slope <- map_dfr(rigal_trends,
~.x %>%
select(siteid, linear_slope),
.id = "response"
)
slope_df <- slope %>%
pivot_wider(names_from = "response", values_from = "linear_slope")
Summary
summary_slope <- slope %>%
group_by(response) %>%
summarise(summ = list(enframe(summary_distribution(linear_slope, na.rm = TRUE)))) %>%
unnest(cols = summ) %>%
pivot_wider(names_from = "name", values_from = "value")
summary_slope %>%
mutate(response = get_var_replacement()[response]) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "bordered", "hover")) %>%
scroll_box(width = "100%", height = "400px")
|
response
|
min
|
1st_quart
|
median
|
2nd_quart
|
max
|
mean
|
sd
|
n
|
n_na
|
frac_na
|
|
Appearance
|
-0.0559441
|
-0.0018271
|
0.0009691
|
0.0065140
|
0.0600000
|
0.0027393
|
0.0098864
|
2678
|
0
|
0
|
|
Chao (binary, similarity)
|
-0.0969697
|
-0.0061754
|
0.0000000
|
0.0000000
|
0.1210526
|
-0.0029650
|
0.0138050
|
2678
|
0
|
0
|
|
Chao evenness
|
-1.0762433
|
-0.0197209
|
0.0000000
|
0.0214427
|
0.8698342
|
0.0000118
|
0.0687473
|
2678
|
0
|
0
|
|
Chao species richness
|
-1.0804000
|
-0.0363727
|
0.0002926
|
0.0473422
|
0.9908720
|
0.0073944
|
0.1387268
|
2678
|
0
|
0
|
|
Chao shannon
|
-0.7782364
|
-0.0270535
|
0.0007389
|
0.0362457
|
0.7822804
|
0.0043759
|
0.0822483
|
2678
|
0
|
0
|
|
Chao simpson
|
-0.6747701
|
-0.0241523
|
0.0007189
|
0.0304959
|
0.5440252
|
0.0028965
|
0.0691378
|
2678
|
0
|
0
|
|
Disappearance
|
-0.0486529
|
-0.0005100
|
0.0000000
|
0.0042135
|
0.0557576
|
0.0017718
|
0.0071868
|
2678
|
0
|
0
|
|
Evenness
|
-0.0830806
|
-0.0060849
|
0.0000000
|
0.0069975
|
0.0817728
|
0.0006221
|
0.0144629
|
2678
|
0
|
0
|
|
SER_a (rel abundance)
|
-0.1059755
|
-0.0146337
|
-0.0032189
|
0.0007927
|
0.1169447
|
-0.0072135
|
0.0183244
|
2678
|
0
|
0
|
|
Horn (binary, similarity)
|
-0.0764433
|
-0.0091583
|
-0.0027316
|
0.0012988
|
0.0859619
|
-0.0041715
|
0.0113065
|
2678
|
0
|
0
|
|
Jaccard (binary, similarity)
|
-0.0781734
|
-0.0115266
|
-0.0037510
|
0.0014597
|
0.0672796
|
-0.0051970
|
0.0135836
|
2678
|
0
|
0
|
|
Jaccard (binary, dissimilarity)
|
-0.0672796
|
-0.0014597
|
0.0037510
|
0.0115266
|
0.0781734
|
0.0051970
|
0.0135836
|
2678
|
0
|
0
|
|
Log species richness
|
-20.6068500
|
-0.9402512
|
0.0498002
|
1.4080463
|
26.9911617
|
0.2947554
|
2.6535099
|
2678
|
0
|
0
|
|
Log total abundance
|
-31.2811462
|
-2.8068822
|
0.4993106
|
4.1780859
|
55.8343888
|
0.8272072
|
6.6756922
|
2678
|
0
|
0
|
|
Nestedness (jaccard)
|
-0.0830556
|
-0.0041642
|
0.0002881
|
0.0074839
|
0.1156349
|
0.0021108
|
0.0137269
|
2678
|
0
|
0
|
|
Shannon
|
-0.2007711
|
-0.0088851
|
0.0004650
|
0.0121768
|
0.1574411
|
0.0016449
|
0.0221596
|
2678
|
0
|
0
|
|
Simpson
|
-0.0829149
|
-0.0045179
|
0.0001071
|
0.0058510
|
0.0800458
|
0.0007642
|
0.0111222
|
2678
|
0
|
0
|
|
Species richness
|
-1.6727273
|
-0.0328546
|
0.0020269
|
0.0536209
|
1.5727273
|
0.0148247
|
0.1422215
|
2678
|
0
|
0
|
|
Total turnover (codyn)
|
-0.0559441
|
-0.0013253
|
0.0032257
|
0.0097964
|
0.0621302
|
0.0045111
|
0.0115567
|
2678
|
0
|
0
|
|
Total abundance
|
-755.8336538
|
-1.8770744
|
0.1129827
|
2.2456287
|
398.4499936
|
1.7786740
|
31.5332673
|
2678
|
0
|
0
|
|
Turnover (jaccard)
|
-0.1223776
|
0.0000000
|
0.0000000
|
0.0077107
|
0.0989899
|
0.0030863
|
0.0134089
|
2678
|
0
|
0
|
bs <- map(var_temporal_trends, ~summary_distribution(slope_df[[.x]]))
names(bs) <- var_temporal_trends
- Jaccard trends: -0.005 (0.014 s.d.)
(-0.01 in @dornelas_assemblage_2014).
- meaning -0.05 % change per decade in average
- Jaccard decreased in 65%
of the sites (79% in @dornelas_assemblage_2014)
- Richness trends:
- +0.3% richness change per year
(2.7 s.d.),
i.e. 3% per decade.
Check trends
p <- map2(rigal_trends, names(rigal_trends),
~.x %>%
ggplot(aes(x = linear_slope)) +
geom_histogram() +
labs(title = .y)
)
plot_grid(plotlist = p)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

na_siteid_slope <- unique(slope[is.na(slope$linear_slope), ]$siteid)
sd_num <- analysis_dataset %>%
group_by(siteid) %>%
summarise(across(where(is.numeric), ~(sd(.x))))
sd_num %>%
filter(siteid %in% na_siteid_slope) %>%
summary(.)
#> siteid species_nb total_abundance log_total_abundance
#> Length:0 Min. : NA Min. : NA Min. : NA
#> Class :character 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Mode :character Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA
#> log_species_nb year latitude longitude main_bas
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> chao_richness chao_shannon chao_simpson chao_evenness hillebrand
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> total appearance disappearance jaccard horn
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> chao total_abundance_int shannon simpson inv_simpson
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> evenness jaccard_dis nestedness turnover first_year
#> Min. : NA Min. : NA Min. : NA Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA Median : NA Median : NA Median : NA
#> Mean :NaN Mean :NaN Mean :NaN Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA Max. : NA Max. : NA Max. : NA
#> last_year span
#> Min. : NA Min. : NA
#> 1st Qu.: NA 1st Qu.: NA
#> Median : NA Median : NA
#> Mean :NaN Mean :NaN
#> 3rd Qu.: NA 3rd Qu.: NA
#> Max. : NA Max. : NA
- Tremendous trends in abundance:
high_abundance_slope <- slope_df %>%
arrange(desc(abs(total_abundance)))
high_abundance_slope
#> # A tibble: 2,678 × 22
#> siteid total_abundance log_total_abundance species_nb log_species_nb
#> <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 S11735 -756. -0.971 0.245 1.31
#> 2 S6466 -472. -7.44 0.309 1.38
#> 3 S5180 398. 3.91 0.207 0.822
#> 4 S5812 -368. 1.23 0.393 1.56
#> 5 S3863 330. 12.6 0.206 1.53
#> 6 S1787 277. 29.3 0.337 2.69
#> 7 S8577 256. 4.78 0.168 0.800
#> 8 S4501 255. 29.4 0.582 10.8
#> 9 S8158 -249. -6.37 0.108 0.562
#> 10 S6074 245. 2.37 0.198 0.908
#> # … with 2,668 more rows, and 17 more variables: chao_richness <dbl>,
#> # chao_shannon <dbl>, chao_simpson <dbl>, chao_evenness <dbl>, jaccard <dbl>,
#> # horn <dbl>, chao <dbl>, hillebrand <dbl>, total <dbl>, appearance <dbl>,
#> # disappearance <dbl>, evenness <dbl>, shannon <dbl>, simpson <dbl>,
#> # jaccard_dis <dbl>, nestedness <dbl>, turnover <dbl>
site_desc_slope_total_abundance <-
high_abundance_slope %>%
.[["siteid"]]
- The most tremendous temporal trends in abundance are located in
Illinois and France
ti <- filtered_dataset$location %>%
filter(siteid %in% site_desc_slope_total_abundance[1:20]) %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 2154)
sp::plot(ne_countries())
sp::plot(st_geometry(ti), add = TRUE, col = "red", pch = 19)

filtered_dataset$abun_rich_op %>%
filter(siteid %in% site_desc_slope_total_abundance[1:20]) %>%
group_by(siteid) %>%
summarise(across(c("total_abundance", "species_nb"), list(median = median, mean = mean)))
#> # A tibble: 20 × 5
#> siteid total_abundance_me… total_abundance_… species_nb_medi… species_nb_mean
#> * <chr> <dbl> <dbl> <dbl> <dbl>
#> 1 S11735 5948 17416. 18 18.7
#> 2 S1507 86 698. 2 10.4
#> 3 S1787 2790. 2498. 16.5 15.7
#> 4 S2595 882 1306. 15 15.9
#> 5 S3445 1791 2112. 5.5 5.31
#> 6 S3863 1827 3353. 13 13.6
#> 7 S3883 755 1842. 10 9.46
#> 8 S4325 644. 899. 5.5 5.58
#> 9 S4501 2111 1924. 16 14.9
#> 10 S4755 2517 2335. 20 20.1
#> 11 S5180 4335 8431. 23 23.3
#> 12 S5392 1416 4011. 19.5 18.9
#> 13 S5812 4968 12216. 24.5 25.2
#> 14 S5987 1470 1630. 18 14.2
#> 15 S6074 5908. 9804. 23.5 23.1
#> 16 S6466 3166 6465. 27 25.7
#> 17 S8158 3610 5140. 24.5 24.7
#> 18 S8577 4474 7435. 28 26.4
#> 19 S8852 4126 5016. 25 22.8
#> 20 S9031 3078 11660. 22 22.2
p_high_slope_abun <- map(site_desc_slope_total_abundance[1:20],
~plot_community_data(
analysis_dataset %>%
filter(siteid == .x) %>%
select(siteid, year, total_abundance),
y = "total_abundance", x = "year",
smoothing_method = "lm"))
plot_grid(plotlist = p_high_slope_abun, ncol = 3)

s11735_com <- analysis_dataset %>%
filter(siteid == site_desc_slope_total_abundance[1]) %>%
arrange(year) %>%
select(siteid, year, total_abundance)
s11735_com
#> # A tibble: 25 × 3
#> siteid year total_abundance
#> <chr> <dbl> <dbl>
#> 1 S11735 1993 66
#> 2 S11735 1994 1366
#> 3 S11735 1995 109088
#> 4 S11735 1996 9997
#> 5 S11735 1997 512
#> 6 S11735 1998 5948
#> 7 S11735 1999 85163
#> 8 S11735 2000 10785
#> 9 S11735 2001 21720
#> 10 S11735 2002 11719
#> # … with 15 more rows
plot_community_data(s11735_com, y = "total_abundance", x = "year",
smoothing_method = "lm")

mean(s11735_com$total_abundance) - sd(s11735_com$total_abundance)
#> [1] -10662.47
mean(s11735_com$total_abundance) + 5 * sd(s11735_com$total_abundance)
#> [1] 157805.5
s11735_pop <- filtered_dataset$measurement %>%
filter(siteid == site_desc_slope_total_abundance[1]) %>%
arrange(year) %>%
select(siteid, year, species, abundance)
s11735_pop %>%
group_by(species) %>%
summarise(abun = sum(abundance)) %>%
arrange(desc(abun))
#> # A tibble: 52 × 2
#> species abun
#> <chr> <dbl>
#> 1 Notropis atherinoides 418706
#> 2 Notropis volucellus 8400
#> 3 Cyprinella spiloptera 1949
#> 4 Lepomis macrochirus 1659
#> 5 Dorosoma cepedianum 651
#> 6 Aplodinotus grunniens 614
#> 7 Morone chrysops 540
#> 8 Pimephales vigilax 527
#> 9 Cyprinus carpio 427
#> 10 Macrhybopsis aestivalis 386
#> # … with 42 more rows
plot_temporal_population(s11735_pop) +
theme(legend.position = "bottom")
#> Warning: Removed 52 rows containing missing values (position_stack).

- Second highest temporal trends in abundance:
s6466_com <- analysis_dataset %>%
filter(siteid == site_desc_slope_total_abundance[2]) %>%
arrange(year) %>%
select(siteid, year, total_abundance)
s6466_com
#> # A tibble: 11 × 3
#> siteid year total_abundance
#> <chr> <dbl> <dbl>
#> 1 S6466 1981 7844
#> 2 S6466 1983 11938
#> 3 S6466 1987 9521
#> 4 S6466 1988 23997
#> 5 S6466 1989 1558
#> 6 S6466 1990 1674
#> 7 S6466 1991 3166
#> 8 S6466 1992 862
#> 9 S6466 1997 2687
#> 10 S6466 1999 1558
#> 11 S6466 2000 6314
plot_community_data(s6466_com, y = "total_abundance", x = "year",
smoothing_method = "lm")

s6466_pop <- filtered_dataset$measurement %>%
filter(siteid == site_desc_slope_total_abundance[2]) %>%
arrange(year) %>%
select(siteid, year, species, abundance)
s6466_pop %>%
group_by(species) %>%
summarise(abun = sum(abundance)) %>%
arrange(desc(abun))
#> # A tibble: 55 × 2
#> species abun
#> <chr> <dbl>
#> 1 Dorosoma cepedianum 43869
#> 2 Aplodinotus grunniens 8914
#> 3 Notropis wickliffi 8637
#> 4 Alosa chrysochloris 1812
#> 5 Morone chrysops 1538
#> 6 Cyprinus carpio 1272
#> 7 Ictalurus punctatus 1192
#> 8 Notropis atherinoides 745
#> 9 Ictiobus bubalus 559
#> 10 Sander canadensis 465
#> # … with 45 more rows
plot_temporal_population(s6466_pop) +
theme(legend.position = "bottom")
#> Warning: Removed 495 rows containing missing values (position_stack).

Relationships among temporal trends
Globally
cor_slope <- slope_df %>%
select(-siteid) %>%
cor(., method = "spearman")
cor_slope
#> total_abundance log_total_abundance species_nb
#> total_abundance 1.000000000 0.830317697 0.277156569
#> log_total_abundance 0.830317697 1.000000000 0.335207093
#> species_nb 0.277156569 0.335207093 1.000000000
#> log_species_nb 0.250125492 0.341531037 0.938150124
#> chao_richness 0.016699737 0.078327544 0.815025195
#> chao_shannon -0.060971568 -0.004315524 0.598119627
#> chao_simpson -0.068607673 -0.030703462 0.486160737
#> chao_evenness -0.046333174 0.001842386 0.158532100
#> jaccard 0.029465496 0.058720565 -0.023342671
#> horn 0.032107656 0.065162610 -0.008994887
#> chao 0.049896277 0.074583303 0.068878719
#> hillebrand 0.049560233 0.095796198 -0.025148726
#> total 0.020626140 0.011421991 0.205181084
#> appearance 0.161650857 0.192922044 0.687912675
#> disappearance -0.161078700 -0.222194738 -0.555815437
#> evenness -0.135468004 -0.106903064 0.191546294
#> shannon -0.030652568 0.031284774 0.595797755
#> simpson -0.036762181 0.012861550 0.475730622
#> jaccard_dis -0.028841716 -0.058365304 0.023323286
#> nestedness 0.004895141 -0.007242935 0.105858491
#> turnover -0.022365349 -0.038795778 -0.052263750
#> log_species_nb chao_richness chao_shannon chao_simpson
#> total_abundance 0.250125492 0.01669974 -0.060971568 -0.06860767
#> log_total_abundance 0.341531037 0.07832754 -0.004315524 -0.03070346
#> species_nb 0.938150124 0.81502519 0.598119627 0.48616074
#> log_species_nb 1.000000000 0.79177426 0.652142713 0.55798584
#> chao_richness 0.791774260 1.00000000 0.715598168 0.58506236
#> chao_shannon 0.652142713 0.71559817 1.000000000 0.94429606
#> chao_simpson 0.557985838 0.58506236 0.944296056 1.00000000
#> chao_evenness 0.238494063 0.17861100 0.627245251 0.67332684
#> jaccard -0.022860296 -0.03144932 -0.016411777 -0.02714028
#> horn -0.002333034 -0.01954284 -0.001651297 -0.01305221
#> chao 0.078669559 0.05198781 0.064378306 0.04836406
#> hillebrand -0.022070175 -0.03353840 -0.068107966 -0.08324337
#> total 0.197786137 0.18733040 0.126450573 0.11238128
#> appearance 0.691348158 0.58925101 0.453612568 0.38451499
#> disappearance -0.577843872 -0.45943156 -0.378710413 -0.30870107
#> evenness 0.289979041 0.32190903 0.722795822 0.77814032
#> shannon 0.652755093 0.67989789 0.944628419 0.91122873
#> simpson 0.555069101 0.55478316 0.873409696 0.89103134
#> jaccard_dis 0.022736545 0.03132624 0.016070172 0.02665117
#> nestedness 0.097239453 0.10838694 0.094533541 0.07288059
#> turnover -0.048541424 -0.05721969 -0.064561731 -0.03107103
#> chao_evenness jaccard horn chao
#> total_abundance -0.046333174 0.029465496 0.032107656 0.04989628
#> log_total_abundance 0.001842386 0.058720565 0.065162610 0.07458330
#> species_nb 0.158532100 -0.023342671 -0.008994887 0.06887872
#> log_species_nb 0.238494063 -0.022860296 -0.002333034 0.07866956
#> chao_richness 0.178610995 -0.031449324 -0.019542841 0.05198781
#> chao_shannon 0.627245251 -0.016411777 -0.001651297 0.06437831
#> chao_simpson 0.673326844 -0.027140280 -0.013052206 0.04836406
#> chao_evenness 1.000000000 0.063506484 0.074803390 0.06164999
#> jaccard 0.063506484 1.000000000 0.989134887 0.69772222
#> horn 0.074803390 0.989134887 1.000000000 0.73627697
#> chao 0.061649988 0.697722217 0.736276969 1.00000000
#> hillebrand 0.014641849 0.332833278 0.347381250 0.24936211
#> total -0.042854829 -0.933845455 -0.920529480 -0.59959876
#> appearance 0.090571956 -0.547707104 -0.527710202 -0.28291096
#> disappearance -0.143240688 -0.642940230 -0.662605243 -0.52031520
#> evenness 0.712515794 0.047910610 0.058183834 0.06970014
#> shannon 0.565547657 -0.005489195 0.009607005 0.07637940
#> simpson 0.595960699 -0.002340240 0.012570159 0.06921570
#> jaccard_dis -0.064111974 -0.999857414 -0.989090989 -0.69677143
#> nestedness -0.008434623 -0.405843668 -0.355504999 -0.33184081
#> turnover -0.037696468 -0.492788826 -0.522609730 -0.27616912
#> hillebrand total appearance disappearance
#> total_abundance 0.04956023 0.02062614 0.16165086 -0.16107870
#> log_total_abundance 0.09579620 0.01142199 0.19292204 -0.22219474
#> species_nb -0.02514873 0.20518108 0.68791268 -0.55581544
#> log_species_nb -0.02207017 0.19778614 0.69134816 -0.57784387
#> chao_richness -0.03353840 0.18733040 0.58925101 -0.45943156
#> chao_shannon -0.06810797 0.12645057 0.45361257 -0.37871041
#> chao_simpson -0.08324337 0.11238128 0.38451499 -0.30870107
#> chao_evenness 0.01464185 -0.04285483 0.09057196 -0.14324069
#> jaccard 0.33283328 -0.93384545 -0.54770710 -0.64294023
#> horn 0.34738125 -0.92052948 -0.52771020 -0.66260524
#> chao 0.24936211 -0.59959876 -0.28291096 -0.52031520
#> hillebrand 1.00000000 -0.31748997 -0.21212267 -0.20672017
#> total -0.31748997 1.00000000 0.72528811 0.48889404
#> appearance -0.21212267 0.72528811 1.00000000 -0.11647853
#> disappearance -0.20672017 0.48889404 -0.11647853 1.00000000
#> evenness -0.05321078 -0.01433782 0.13964393 -0.17573304
#> shannon -0.07563072 0.11786971 0.44713913 -0.38733767
#> simpson -0.07581108 0.08965621 0.36745811 -0.32242970
#> jaccard_dis -0.33261575 0.93400650 0.54788050 0.64308627
#> nestedness -0.07413472 0.38337388 0.24465363 0.08227401
#> turnover -0.21287472 0.47223951 0.30131921 0.50349434
#> evenness shannon simpson jaccard_dis
#> total_abundance -0.135468004 -0.030652568 -0.036762181 -0.028841716
#> log_total_abundance -0.106903064 0.031284774 0.012861550 -0.058365304
#> species_nb 0.191546294 0.595797755 0.475730622 0.023323286
#> log_species_nb 0.289979041 0.652755093 0.555069101 0.022736545
#> chao_richness 0.321909029 0.679897895 0.554783164 0.031326242
#> chao_shannon 0.722795822 0.944628419 0.873409696 0.016070172
#> chao_simpson 0.778140322 0.911228727 0.891031340 0.026651172
#> chao_evenness 0.712515794 0.565547657 0.595960699 -0.064111974
#> jaccard 0.047910610 -0.005489195 -0.002340240 -0.999857414
#> horn 0.058183834 0.009607005 0.012570159 -0.989090989
#> chao 0.069700135 0.076379395 0.069215704 -0.696771428
#> hillebrand -0.053210777 -0.075630723 -0.075811084 -0.332615747
#> total -0.014337815 0.117869711 0.089656212 0.934006505
#> appearance 0.139643932 0.447139132 0.367458109 0.547880499
#> disappearance -0.175733040 -0.387337667 -0.322429697 0.643086266
#> evenness 1.000000000 0.784047749 0.866977324 -0.048397766
#> shannon 0.784047749 1.000000000 0.962316417 0.005099809
#> simpson 0.866977324 0.962316417 1.000000000 0.001943911
#> jaccard_dis -0.048397766 0.005099809 0.001943911 1.000000000
#> nestedness 0.009438247 0.090073607 0.065840971 0.405983665
#> turnover -0.053642052 -0.071596893 -0.052102504 0.492799673
#> nestedness turnover
#> total_abundance 0.004895141 -0.02236535
#> log_total_abundance -0.007242935 -0.03879578
#> species_nb 0.105858491 -0.05226375
#> log_species_nb 0.097239453 -0.04854142
#> chao_richness 0.108386938 -0.05721969
#> chao_shannon 0.094533541 -0.06456173
#> chao_simpson 0.072880587 -0.03107103
#> chao_evenness -0.008434623 -0.03769647
#> jaccard -0.405843668 -0.49278883
#> horn -0.355504999 -0.52260973
#> chao -0.331840814 -0.27616912
#> hillebrand -0.074134723 -0.21287472
#> total 0.383373875 0.47223951
#> appearance 0.244653628 0.30131921
#> disappearance 0.082274006 0.50349434
#> evenness 0.009438247 -0.05364205
#> shannon 0.090073607 -0.07159689
#> simpson 0.065840971 -0.05210250
#> jaccard_dis 0.405983665 0.49279967
#> nestedness 1.000000000 -0.47739313
#> turnover -0.477393127 1.00000000
png(here("doc", "fig", "p_cor_slope.png"))
corrplot::corrplot(cor_slope[c("chao_richness", "species_nb", "log_species_nb",
"total_abundance", "log_total_abundance"), ], type = "upper", tl.srt = 45, diag = FALSE)
dev.off()
#> png
#> 2
png(here("doc", "fig", "p_cor_slope_tot.png"))
corrplot::corrplot(cor_slope, type = "upper", tl.srt = 35, diag = FALSE)
dev.off()
#> png
#> 2
corrplot::corrplot(cor_slope[c("chao_richness", "species_nb", "log_species_nb",
"total_abundance", "log_total_abundance"), ], type = "upper", tl.srt = 45, diag = FALSE)

corrplot::corrplot(cor_slope, type = "upper", tl.srt = 35, diag = FALSE)

ti <- expand.grid(
resp1 = unique(slope$response),
resp2 = unique(slope$response)
) %>%
filter(resp2 != resp1) %>%
filter(
resp1 %in% c("chao_richness", "species_nb", "log_species_nb",
"total_abundance", "log_total_abundance")) %>%
mutate_all(as.character) %>%
arrange(resp1)
test <- map2(ti$resp1, ti$resp2,
function(x, y) {
bi <- slope %>%
filter(response %in% c(x, y)) %>%
select(siteid, response, linear_slope) %>%
pivot_wider(names_from = "response", values_from = "linear_slope")
return(bi)
}
)
p_trends_trends <- map(test, function(x) {
l <- colnames(x)
x %>%
ggplot(aes(x = !!sym(l[2]), y = !!sym(l[3]))) +
geom_point() +
geom_smooth(method = "loess") +
labs(
x = get_var_replacement()[l[2]],
y = get_var_replacement()[l[3]]
)
})
names(p_trends_trends) <- map_chr(test, ~colnames(.x)[2])
Total abundance
plot_grid(
plotlist = p_trends_trends[names(p_trends_trends) %in% "total_abundance"],
ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

Species richness
plot_grid(
plotlist = p_trends_trends[names(p_trends_trends) %in% "species_nb"],
ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

Log species richness
plot_grid(
plotlist = p_trends_trends[names(p_trends_trends) %in% "log_species_nb"],
ncol = 3
)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

- Trends of total abundance are linked to the one of species richness but not to
shannon and simpson, suggesting that the trends of abundance are more link to
trends of rare species. Logic: gain and loss of species should concern rare
species.
- Species richness decrease linked to disappearance, but not increase
- Species richness increase linked to appearance, but not decrease
Turnover vs nestedness:
test <- slope_df %>%
mutate(
nes_sup_tu = nestedness > turnover,
nes_inf_tu = nestedness < turnover,
rich_inc = species_nb >= 0
)
sum(test$nes_sup_tu)
#> [1] 1224
sum(test$nes_inf_tu)
#> [1] 1318
test %>%
tabyl(rich_inc, nes_sup_tu) %>%
adorn_totals("row") %>%
adorn_title(placement = "top")
#> nes_sup_tu
#> rich_inc FALSE TRUE
#> FALSE 733 499
#> TRUE 721 725
#> Total 1454 1224
median_site_rich <- analysis_dataset %>%
select(c("siteid", starts_with("chao_"), "total_abundance")) %>%
group_by(siteid) %>%
summarise(across(where(is.double), median)) %>%
rename_with(~paste0("med_", .x), where(is.double))
slope_df_med <- slope_df %>%
left_join(median_site_rich, by = "siteid") %>%
mutate(
nes_sup_tu = nestedness > turnover,
nes_inf_tu = nestedness < turnover,
rich_inc = species_nb >= 0
)
slope_df_med %>%
ggplot(aes(
x = log_species_nb,
y = jaccard_dis,
shape = nes_sup_tu,
color = nes_sup_tu
)) +
geom_point() +
geom_smooth(method = "gam")
#> `geom_smooth()` using formula 'y ~ s(x, bs = "cs")'

p_baselga_jaccard <- map(c("nestedness", "turnover"),
~slope_df_med %>%
ggplot(aes(
x = jaccard_dis,
y = !!sym(.x),
color = log(med_chao_richness))) +
geom_point() +
scale_color_viridis() +
geom_abline(intercept = 0, slope = 1) +
geom_smooth(method = "loess")
)
p_baselga_jaccard[[3]] <-
slope_df_med %>%
ggplot(aes(
x = jaccard_dis,
y = nestedness - turnover,
color = log(med_chao_richness))) +
geom_point() +
scale_color_viridis() +
geom_abline(intercept = 0, slope = 1) +
geom_smooth(method = "loess")
plot_grid(plotlist = p_baselga_jaccard, ncol = 3)
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'
#> `geom_smooth()` using formula 'y ~ x'

slope_df_med %>%
ggplot(aes(x = med_chao_richness, y = jaccard_dis)) +
geom_point()

Look at strong turnover trends
high_jaccard_slope <-
slope_df %>%
select(siteid, total_abundance, species_nb, jaccard) %>%
arrange(desc(jaccard))
high_jaccard_slope
#> # A tibble: 2,678 × 4
#> siteid total_abundance species_nb jaccard
#> <chr> <dbl> <dbl> <dbl>
#> 1 S450 120. 0.870 0.0673
#> 2 S6544 0.257 -0.0769 0.0559
#> 3 S9623 -3.77 -0.0266 0.0556
#> 4 S3217 0.0387 -0.0473 0.0496
#> 5 S6550 1.92 0.116 0.0490
#> 6 S4875 28.3 0.145 0.0461
#> 7 S2831 6.25 -0.0455 0.0420
#> 8 S9060 -0.817 -0.0303 0.0417
#> 9 S5514 0.690 -0.0576 0.0403
#> 10 S2742 1.08 -0.00909 0.0394
#> # … with 2,668 more rows
site_desc_jaccard_slope <-
high_jaccard_slope %>%
.[["siteid"]]
ti <- filtered_dataset$location %>%
filter(siteid %in% site_desc_jaccard_slope[1:20]) %>%
st_as_sf(coords = c("longitude", "latitude"), crs = 2154)
sp::plot(ne_countries())
sp::plot(st_geometry(ti), add = TRUE, col = "red", pch = 19)

filtered_dataset$abun_rich_op %>%
filter(siteid %in% site_desc_jaccard_slope[1:20]) %>%
group_by(siteid) %>%
summarise(across(c("total_abundance", "species_nb"),
list(median = median, mean = mean))) %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "bordered", "hover")) %>%
scroll_box(width = "100%", height = "400px")
|
siteid
|
total_abundance_median
|
total_abundance_mean
|
species_nb_median
|
species_nb_mean
|
|
S11001
|
29.5700
|
33.575000
|
4.0
|
4.181818
|
|
S1574
|
323.0000
|
382.245454
|
4.0
|
3.272727
|
|
S2430
|
40.6520
|
39.779385
|
1.0
|
1.461539
|
|
S2742
|
41.9000
|
57.090909
|
2.0
|
2.000000
|
|
S2831
|
170.5000
|
180.154546
|
1.0
|
1.272727
|
|
S3044
|
40.5000
|
50.280000
|
1.0
|
1.100000
|
|
S3217
|
7.7000
|
8.793333
|
3.0
|
2.800000
|
|
S450
|
877.5000
|
744.666667
|
12.5
|
9.166667
|
|
S4875
|
290.5000
|
258.900000
|
13.0
|
11.200000
|
|
S5046
|
546.5000
|
596.100000
|
15.5
|
11.200000
|
|
S5514
|
25.9000
|
23.580000
|
1.5
|
1.700000
|
|
S560
|
44.1500
|
45.183333
|
3.0
|
3.416667
|
|
S6544
|
15.0000
|
18.061539
|
1.0
|
1.307692
|
|
S6550
|
6.6500
|
14.557143
|
2.0
|
2.214286
|
|
S6783
|
10.1000
|
12.246667
|
1.0
|
1.200000
|
|
S6936
|
22.8000
|
26.625000
|
1.0
|
1.333333
|
|
S7806
|
10.2000
|
10.950000
|
1.0
|
1.222222
|
|
S8447
|
10.7000
|
11.000000
|
4.0
|
3.750000
|
|
S9060
|
21.8185
|
21.676300
|
1.5
|
1.500000
|
|
S9623
|
61.0370
|
54.827200
|
2.5
|
2.500000
|
p_high_jaccard_slope <- map(site_desc_jaccard_slope[1:20],
~plot_community_data(
analysis_dataset %>%
filter(siteid == .x) %>%
select(siteid, year, jaccard),
y = "jaccard", x = "year",
smoothing_method = "lm"))
plot_grid(plotlist = p_high_jaccard_slope, ncol = 3)
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).

tar_load(vegdist_turnover_c)
- Hum, do not understand why community are more similar at the end of the time-steps
p_pop_high_jaccard <-
map(site_desc_jaccard_slope[1:20],
~plot_temporal_population(
filtered_dataset$measurement %>%
filter(siteid == .x)
) +
theme(legend.position = "none")
)
plot_grid(plotlist = p_pop_high_jaccard, ncol = 3)
#> Warning: Removed 24 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 77 rows containing missing values (position_stack).
#> Warning: Removed 8 rows containing missing values (position_stack).
#> Warning: Removed 36 rows containing missing values (position_stack).
#> Warning: Removed 4 rows containing missing values (position_stack).
#> Warning: Removed 3 rows containing missing values (position_stack).
#> Warning: Removed 8 rows containing missing values (position_stack).
#> Warning: Removed 60 rows containing missing values (position_stack).
#> Warning: Removed 2 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 26 rows containing missing values (position_stack).

- Lets see communities that becomes more dissimilar:
p_jaccard_low_jaccard <-
map(site_desc_jaccard_slope[(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
~plot_community_data(
analysis_dataset %>%
filter(siteid == .x) %>%
select(siteid, year, jaccard),
y = "jaccard", x = "year",
smoothing_method = "lm"))
plot_grid(plotlist = p_jaccard_low_jaccard, ncol = 3)
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).
#> Warning: Removed 1 rows containing non-finite values (stat_smooth).
#> Warning: Removed 1 rows containing missing values (geom_point).

p_pop_low_jaccard <-
map(site_desc_jaccard_slope[
(length(site_desc_jaccard_slope)-20):length(site_desc_jaccard_slope)],
~plot_temporal_population(
filtered_dataset$measurement %>%
filter(siteid == .x)
) +
theme(legend.position = "none")
)
plot_grid(plotlist = p_pop_low_jaccard, ncol = 3)
#> Warning: Removed 60 rows containing missing values (position_stack).
#> Warning: Removed 24 rows containing missing values (position_stack).
#> Warning: Removed 312 rows containing missing values (position_stack).
#> Warning: Removed 4 rows containing missing values (position_stack).
#> Warning: Removed 9 rows containing missing values (position_stack).
#> Warning: Removed 30 rows containing missing values (position_stack).
#> Warning: Removed 63 rows containing missing values (position_stack).
#> Warning: Removed 28 rows containing missing values (position_stack).
#> Warning: Removed 4 rows containing missing values (position_stack).
#> Warning: Removed 36 rows containing missing values (position_stack).
#> Warning: Removed 60 rows containing missing values (position_stack).
#> Warning: Removed 4 rows containing missing values (position_stack).
#> Warning: Removed 18 rows containing missing values (position_stack).
#> Warning: Removed 12 rows containing missing values (position_stack).
#> Warning: Removed 27 rows containing missing values (position_stack).
#> Warning: Removed 126 rows containing missing values (position_stack).

Comparison rigal, lm
rigal_trends_df <- map2_dfr(
rigal_trends, names(rigal_trends),
~.x %>% mutate(variable = .y)
) %>%
select(variable, siteid, linear_slope) %>%
pivot_wider(names_from = "variable", values_from = "linear_slope") %>%
mutate(type = "rigal")
slope_comp <- list()
slope_comp$rigal <- rigal_trends_df[order(rigal_trends_df$siteid), ]
slope_comp$simple_lm <- slope_df[order(slope_df$siteid), ]
stopifnot(all(slope_comp$rigal$siteid == slope_comp$simple_lm$siteid))
slope_comp$diff <- tibble(
siteid = slope_comp$rigal$siteid)
for (i in var_temporal_trends) {
slope_comp$diff[[i]] <- slope_comp$rigal[[i]] - slope_comp$simple_lm[[i]]
}
old_par <- par()
par(mfrow = c(3, 2))
for (i in var_temporal_trends) {
plot(
slope_comp$rigal[[i]] ~
slope_comp$simple_lm[[i]],
ylab = "Rigal (polynomial degree 2)",
xlab = "Simple LM"
)
abline(0, 1)
title(i)
}



par(old_par)
#> Warning in par(old_par): le paramètre graphique "cin" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cra" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "csi" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cxy" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "din" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "page" ne peut être changé

# Compare Rigal and simple lm for few sites
plot_rigal_raw_data <- function(
site = NULL,
response = NULL,
dataset = NULL,
rigal_trends = NULL
) {
raw_data <- dataset %>%
filter(siteid == site)
rigal_resp <- rigal_trends[[response]]
rigal_resp_site <- rigal_resp[rigal_resp$siteid == site, ]
plot(raw_data[[response]]~raw_data[["year"]], ylab = response, xlab = "Year")
abline(rigal_resp_site$intercept, rigal_resp_site$linear_slope, col
= "red")
abline(lm(raw_data[[response]]~raw_data[["year"]]), col = "green")
title(main = site)
}
plot_rigal_raw_data(site = "S2672", response = "log_total_abundance", dataset =
analysis_dataset, rigal_trends = rigal_trends)

slope_comp$diff <- slope_comp$diff %>%
arrange(desc(abs(log_total_abundance)))
old_par <- par()
par(mfrow = c(3, 2))
for (i in slope_comp$diff[1:6,]$siteid) {
plot_rigal_raw_data(
site = i,
response = "log_total_abundance",
dataset = analysis_dataset,
rigal_trends = rigal_trends
)
}

par(old_par)
#> Warning in par(old_par): le paramètre graphique "cin" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cra" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "csi" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "cxy" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "din" ne peut être changé
#> Warning in par(old_par): le paramètre graphique "page" ne peut être changé
Classification of temporal trends
classification_station <- map2(rigal_trends, names(rigal_trends),
~.x %>%
ggplot(aes(x = shape_class)) +
geom_bar() +
labs(x = "Trajectory classification", y = "Nb station",
title = get_var_replacement()[.y]) +
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))
)
plot_grid(plotlist = classification_station)

p_classif_abun_rich <- plot_grid(
plotlist = classification_station[
names(classification_station)[str_detect(names(classification_station),
"abundance|richness|species")]])
save_plot(
here("doc", "fig", "classif_trends_abun_rich.png"),
p_classif_abun_rich
)
p_classif_hill_evenness <- plot_grid(
plotlist = classification_station[
names(classification_station)[str_detect(names(classification_station),
"shannon|simpson|evennes")]])
save_plot(
here("doc", "fig", "classif_trends_hill_evenness.png"),
p_classif_hill_evenness
)
p_classif_turnover <- plot_grid(
plotlist = classification_station[
names(classification_station)[names(classification_station) %in%
c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
"disappearance")]])
save_plot(
here("doc", "fig", "classif_trends_classif_hill_evenness.png"),
p_classif_turnover
)
ti <- map_dfr(rigal_trends, ~tabyl_df(x = .x, group = "direction"),
.id = "response"
)
ti %>%
filter(direction != "Total") %>%
select(response, direction, percent) %>%
mutate(response = str_replace_all(response, get_var_replacement())) %>%
pivot_wider(names_from = "direction", values_from = "percent") %>%
kable() %>%
kable_styling(bootstrap_options = c("striped", "bordered", "hover"))
|
response
|
stable
|
increase
|
decrease
|
|
Total abundance
|
75.8%
|
14.3%
|
9.9%
|
|
Log Total Turnover (jaccard) (codyn) abundance
|
73.3%
|
16.5%
|
10.2%
|
|
Species richness
|
77.7%
|
13.2%
|
7.4%
|
|
Log species richness
|
76.2%
|
13.2%
|
6.6%
|
|
Chao species richness
|
79.8%
|
10.6%
|
8.2%
|
|
Chao Shannon
|
76.5%
|
12.1%
|
7.7%
|
|
Chao Simpson
|
78.4%
|
10.5%
|
7.4%
|
|
Chao Evenness
|
82.8%
|
7.5%
|
5.9%
|
|
Jaccard (binary, similarity)
|
76.8%
|
2.8%
|
17.5%
|
|
Horn (binary, similarity)
|
77.3%
|
2.7%
|
17.2%
|
|
Chao (binary, similarity)
|
71.5%
|
2.3%
|
8.3%
|
|
SER_a (rel abundance)
|
75.3%
|
4.1%
|
18.6%
|
|
Total Turnover (jaccard) (codyn)
|
73.7%
|
18.3%
|
3.0%
|
|
Appearance
|
70.2%
|
15.6%
|
3.8%
|
|
disAppearance
|
64.0%
|
9.5%
|
3.0%
|
|
Evenness
|
78.8%
|
9.4%
|
8.1%
|
|
Shannon
|
75.7%
|
12.2%
|
8.3%
|
|
Simpson
|
77.9%
|
10.8%
|
7.7%
|
|
Jaccard (binary, dissimilarity)
|
74.7%
|
17.4%
|
2.8%
|
|
Nestedness (jaccard)
|
80.3%
|
9.3%
|
5.3%
|
|
Turnover (jaccard)
|
57.4%
|
8.3%
|
1.7%
|
#3.3 Log-linear model:logYi=α+βXi+εiIn
#the log-linear model, the literal interpretation of the estimated
#coefficientˆβis that a one-unitincrease inXwill produce an expected increase in
#logYofˆβunits.
#In terms ofYitself, this
#meansthat the expected value ofYis multiplied byeˆβ. So in terms of
#effects of changes inXonY(unlogged):•Each 1-unit increase inXmultiplies the
#expected value ofYbyeˆβ.•To compute the effects onYof another change inXthan an
#increase of one unit, call thischangec, we need to includecin the exponent. The
#effect of ac-unit increase inXis tomultiply the expected value ofYbyecˆβ. So the
#effect for a 5-unit increase inXwould bee5ˆβ.•For small values ofˆβ,
#approximatelyeˆβ≈1+ˆβ. We can use this for the following approxima-tion for a
#quick interpretation of the coefficients: 100·ˆβis the expected percentage
#changeinYfor a unit increase inX. For instance forˆβ=.06,e.06≈1.06, so a 1-unit
#change inXcorresponds to (approximately) an expected increase inYof 6%
#https://data.library.virginia.edu/interpreting-log-transformations-in-a-linear-model/
hist_linear_plot <- map2(rigal_trends, names(rigal_trends),
~.x %>%
ggplot(aes(x = linear_slope)) +
geom_histogram() +
labs(x = "Linear slope", y = "Number of station",
title = get_var_replacement()[.y]) +
ylim(c(0, 1000))
)
plot_grid(plotlist = hist_linear_plot)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).

type_response <- list(
abun_rich = str_detect(names(classification_station),
"abundance|richness|species"),
hill_evenness = str_detect(names(classification_station),
"shannon|simpson|evennes"),
turnover = names(classification_station) %in%
c("jaccard", "horn", "chao", "hillebrand", "total", "appearance",
"disappearance")
)
map2(type_response, names(type_response),
function(x, y) {
p <- plot_grid(
plotlist = hist_linear_plot[
names(hist_linear_plot)[x]
]
)
save_plot(
here("doc", "fig", paste0("hist_linear_slope_", y, ".png")),
p
)
}
)
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#> Warning: Removed 1 rows containing missing values (geom_bar).
#> $abun_rich
#> NULL
#>
#> $hill_evenness
#> NULL
#>
#> $turnover
#> NULL
Spatialization of temporal trends
rigal_trends_df_loc <-
rigal_trends_df %>%
left_join(
select(filtered_dataset$location, siteid, ecoregion, country),
by = c("siteid"))
- No obvious geographic patterns
rigal_trends_df_loc %>%
ggplot(aes(x = country, y = log_species_nb, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc %>%
ggplot(aes(x = country, y = log_total_abundance, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc %>%
ggplot(aes(x = country, y = total, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc %>%
ggplot(aes(x = country, y = jaccard, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc %>%
ggplot(aes(x = country, y = jaccard_dis, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc %>%
ggplot(aes(x = country, y = hillebrand, fill = ecoregion)) +
geom_boxplot()

rigal_trends_df_loc <- filtered_dataset$location %>%
left_join(rigal_trends_df, by = "siteid") %>%
st_as_sf(coords = c("longitude", "latitude"),
crs = 4326)
library(colorspace)
div_pal <- colorspace::sequential_hcl(11, "plasma")
world <- ne_countries(returnclass = "sf")
bb <- st_bbox(rigal_trends_df_loc)
map_world_trends(df = rigal_trends_df_loc, y = "log_species_nb")

map_world_trends(df = rigal_trends_df_loc, y = "jaccard")

map_world_trends(df = rigal_trends_df_loc, y = "log_total_abundance")

rigal_trends_df %>%
summarise(across(where(is.numeric), median, na.rm = T)) %>%
unlist()
#> total_abundance log_total_abundance species_nb log_species_nb
#> 1.129827e-01 4.993106e-01 2.026907e-03 4.980022e-02
#> chao_richness chao_shannon chao_simpson chao_evenness
#> 2.926352e-04 7.389185e-04 7.189205e-04 0.000000e+00
#> jaccard horn chao hillebrand
#> -3.750972e-03 -2.731569e-03 -1.598283e-17 -3.218861e-03
#> total appearance disappearance evenness
#> 3.225717e-03 9.691423e-04 0.000000e+00 0.000000e+00
#> shannon simpson jaccard_dis nestedness
#> 4.650060e-04 1.071122e-04 3.750972e-03 2.880739e-04
#> turnover
#> 0.000000e+00
median(exp(rigal_trends_df$log_total_abundance) - 1) * 100
#> [1] 64.75967
sd(exp(rigal_trends_df$log_total_abundance) - 1) * 100
#> [1] 3.425009e+24
median(exp(na.omit(rigal_trends_df$log_species_nb)) - 1) * 100
#> [1] 5.106305
sd(exp(na.omit(rigal_trends_df$log_species_nb)) - 1) * 100
#> [1] 1.019932e+12
Sensibility analysis of temporal trends
tar_load(c(rigal_trends_avg3y))
names(rigal_trends_avg3y) <- var_temporal_trends
tar_load(c(analysis_dataset_avg3y, measurement_avg3y, measurement))
tar_load(c(turnover_avg3y_c))
tar_load(c(vegdist_turnover_c))
tar_load(c(hillnb_avg3y))
slope_avg3y <- map_dfr(rigal_trends_avg3y[
! names(rigal_trends_avg3y) %in% c("jaccard_dis", "nestedness", "turnover")], ~.x %>%
select(siteid, linear_slope),
.id = "response"
)
- Number of sites that are in the different datasets:
u_site <- map(list(classic = slope, avg_3y = slope_avg3y), ~unique(.x$siteid))
sum(u_site[[1]] %in% u_site[[2]])
#> [1] 2678
sum(u_site[[2]] %in% u_site[[1]])
#> [1] 2678
comparison_baseline <- rbind(
slope %>%
mutate(type = "classic") %>%
filter(siteid %in% u_site[["avg_3y"]]),
slope_avg3y %>%
mutate(type = "avg_3y")
) %>%
pivot_wider(names_from = type, values_from = linear_slope)
- It is not too bad, but we can see some differences:
- slope negatively biaised for
log_species_nb in
3 years baseline
- Bunch of trends in appearance/disappearance trends disappears with
3 years baseline
- But note that in three years baseline, there are three points less data
comparison_baseline %>%
ggplot(aes(x = classic, y = avg_3y)) +
geom_point() +
geom_abline(intercept = 0, slope = 1) +
facet_wrap(~response, scale = "free")
#> Warning: Removed 8034 rows containing missing values (geom_point).

Reproducibility
Reproducibility receipt
## datetime
Sys.time()
#> [1] "2022-01-20 10:22:18 CST"
## repository
if(requireNamespace('git2r', quietly = TRUE)) {
git2r::repository()
} else {
c(
system2("git", args = c("log", "--name-status", "-1"), stdout = TRUE),
system2("git", args = c("remote", "-v"), stdout = TRUE)
)
}
#> Local: main /home/alain/Documents/post-these/isu/RivFishTimeBiodiversityFacets
#> Head: [3695d10] 2022-01-20: fix hillebrand: change entry matrix
## session info
sessionInfo()
#> R version 4.0.5 (2021-03-31)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 10 (buster)
#>
#> Matrix products: default
#> BLAS: /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRblas.so
#> LAPACK: /home/alain/.Renv/versions/4.0.5/lib/R/lib/libRlapack.so
#>
#> locale:
#> [1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
#> [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
#> [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
#>
#> attached base packages:
#> [1] parallel stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] colorspace_2.0-0 betapart_1.5.4 untb_1.7-4
#> [4] iNEXT_2.0.20 factoextra_1.0.7 ade4_1.7-16
#> [7] ecocom_0.0.0.9000 inlabru_2.3.1 INLA_21.11.22
#> [10] sp_1.4-5 foreach_1.5.1 Matrix_1.3-2
#> [13] future_1.21.0 vegan_2.5-7 lattice_0.20-41
#> [16] permute_0.9-5 codyn_2.0.5 janitor_2.1.0
#> [19] viridis_0.5.1 viridisLite_0.3.0 cowplot_1.1.1
#> [22] rnaturalearthdata_0.1.0 rnaturalearth_0.1.0 sf_1.0-4
#> [25] rmarkdown_2.11 scales_1.1.1 kableExtra_1.3.1
#> [28] here_1.0.1 lubridate_1.7.9.2 magrittr_2.0.1
#> [31] forcats_0.5.1 stringr_1.4.0 dplyr_1.0.4
#> [34] purrr_0.3.4 readr_2.1.1 tidyr_1.1.2
#> [37] tibble_3.1.6 ggplot2_3.3.3 tidyverse_1.3.0
#> [40] tarchetypes_0.3.2 targets_0.8.1 conflicted_1.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] snow_0.4-4 readxl_1.3.1 backports_1.2.1
#> [4] fastmatch_1.1-3 corrplot_0.84 plyr_1.8.6
#> [7] igraph_1.2.6 splines_4.0.5 gmp_0.6-2.1
#> [10] listenv_0.8.0 digest_0.6.27 htmltools_0.5.1.1
#> [13] fansi_0.5.0 memoise_2.0.0 cluster_2.1.1
#> [16] tzdb_0.2.0 openxlsx_4.2.3 globals_0.14.0
#> [19] modelr_0.1.8 prettyunits_1.1.1 rvest_0.3.6
#> [22] ggrepel_0.9.1 haven_2.3.1 xfun_0.28
#> [25] rgdal_1.5-28 callr_3.7.0 crayon_1.4.2
#> [28] jsonlite_1.7.2 ape_5.5 iterators_1.0.13
#> [31] glue_1.5.1 gtable_0.3.0 webshot_0.5.2
#> [34] car_3.0-10 abind_1.4-5 emo_0.0.0.9000
#> [37] DBI_1.1.1 rstatix_0.7.0 Rcpp_1.0.6
#> [40] progress_1.2.2 magic_1.5-9 units_0.6-7
#> [43] foreign_0.8-81 httr_1.4.2 RColorBrewer_1.1-2
#> [46] wk_0.5.0 ellipsis_0.3.2 partitions_1.10-4
#> [49] pkgconfig_2.0.3 farver_2.0.3 sass_0.3.1
#> [52] dbplyr_2.1.0 utf8_1.2.2 polynom_1.4-0
#> [55] reshape2_1.4.4 tidyselect_1.1.1 labeling_0.4.2
#> [58] rlang_0.4.12 munsell_0.5.0 cellranger_1.1.0
#> [61] tools_4.0.5 cachem_1.0.4 cli_3.1.0
#> [64] generics_0.1.0 mathjaxr_1.4-0 broom_0.7.4
#> [67] geometry_0.4.5 evaluate_0.14 fastmap_1.1.0
#> [70] yaml_2.2.1 processx_3.5.2 knitr_1.36
#> [73] fs_1.5.1 zip_2.1.1 s2_1.0.7
#> [76] nlme_3.1-152 mime_0.11 xml2_1.3.2
#> [79] compiler_4.0.5 rstudioapi_0.13 curl_4.3.2
#> [82] e1071_1.7-4 ggsignif_0.6.1 reprex_1.0.0
#> [85] bslib_0.2.4 stringi_1.7.6 highr_0.9
#> [88] ps_1.6.0 Brobdingnag_1.2-6 rgeos_0.5-5
#> [91] classInt_0.4-3 vctrs_0.3.8 pillar_1.6.4
#> [94] lifecycle_1.0.1 furrr_0.2.2 jquerylib_0.1.3
#> [97] data.table_1.13.6 R6_2.5.1 bookdown_0.24
#> [100] KernSmooth_2.23-18 gridExtra_2.3 rio_0.5.16
#> [103] parallelly_1.23.0 codetools_0.2-18 rcdd_1.5
#> [106] MASS_7.3-53.1 assertthat_0.2.1 picante_1.8.2
#> [109] rprojroot_2.0.2 withr_2.4.3 doSNOW_1.0.19
#> [112] mgcv_1.8-34 hms_1.1.1 grid_4.0.5
#> [115] class_7.3-18 snakecase_0.11.0 carData_3.0-4
#> [118] git2r_0.29.0 ggpubr_0.4.0 itertools_0.1-3